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Abstract. 

In this paper, we derive consistent shallow water equations for bi-layer flows 
of Newtonian fluids flowing down a ramp. We carry out a complete spectral 
analysis of steady flows in the low frequency regime and show the occurence 
of hydrodynamic instabilities, so called roll-waves, when steady flows are 
unstable. 

1 Introduction 

This paper is devoted to the analysis of the gravity driven motion of a super- 
position of two immiscible Newtonian fluids flowing down an inclined plane. 
Such systems can describe a lot of situations in geophysics and engineering: 
mud flows, submarine avalanches, transport of mass, heat and momentum 
in chemical technology, coating layers in photography. For this latter appli- 
cation, the formation of waves is highly indesirable. It is then an important 
problem to study the stability of such multiple-layer flows. Linear stability 
analysis was addressed in many papers: see e.g. [8] , [9] , [12] . But these studies 
did not provide a model to describe the nonlinear waves. Indeed, modeling 
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such systems is a hard problem both from the mathematical and numerical 
viewpoint: in particular, one has to deal with two free surfaces: one at the 
fluid interface and the other one at the interface between fluid and gas. 

Here, we consider the particular situation where two thin fluids are flowing 
down a ramp. This means that the characteristic depths of the fluids are 
much smaller than the characteristic length of the flow in the downstream 
direction. We take advantage of the thinness of the layers to write a reduced 
system of equations which will contain all the physical ingredients that are 
relevant to describe the dynamics of such flows. A similar strategy was devel- 
oped by Kliakhandler in pTj: a system of Kuramoto Sivashinsky equations 
is derived from the full Navier Stokes equations in the presence of surface 
tension. Using this approach, the author analysed the spectral stability of 
two-layered thin film flows and considered in particular the interaction be- 
tween convection and each relevant physical term: buoyancy, inertia and 
capillarity. In particular, it is proved that in some parameter regime, the 
convection can stabilize an unstable density stratification. Though this spec- 
tral analysis highlights the role of each term (buoyancy, inertia, capillarity), 
it is not complete since one has to consider the interaction between all the 
relevant terms. In particular, the competition between inertia and buoyancy 
is the source of hydrodynamical instabilities in shallow waters. Moreover, 
the derivation of Kuramoto Sivashnisky equations is usually limited to small 
amplitude motions. The purpose of this paper is to obtain a system of shal- 
low water equations which is consistent with Navier-Stokes equations in the 
regime of shallow waters. In order to derive such a system, we follow the 
methodology introduced by Vila [19] and justified rigorously by Bresch and 
Noble |1] for a single fluid layer. As a byproduct, the system of shallow wa- 
ter equations is relevant to study the linear stability of steady flows in the 
low frequency regime. We will complete the spectral analysis of [TT] in some 
particular cases (stable/unstable mass stratification, viscous stratification). 
We prove that the system of Kuramato Sivashinsky equations in [TT] is also 
contained in our model in some specific regimes. Finally, we use shallow 
water equations to describe nonlinear waves when steady flows are spectrally 
unstable. In particular, we show the occurence of well known hydrodynamic 
instabilities, so called roll-waves, which can appear either at fluid interface 
and free surface or only at the fluid interface. 

The paper is organized as follows. In section [21 we describe bi-layer flows 
in the shallow water scaling and compute an expansion of the velocity and 
pressure field in this regime. From this asymptotic analysis, we find the 
system of Kuramoto Sivashinsky equations of [TT] and study the spectral 
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stability of steady flows in the low frequency regime. Then we derive shallow 
water equations. In section [31 we prove the existence of small amplitude 
roll-waves when the steady flow is unstable. Finally, we compute numeri- 
cally large amplitude roll-waves through direct numerical simulations of the 
shallow water equations. 



2 Shallow Water Eqs. for bi- layer flows 

In this section, we show how to expand solutions to Navier-Stokes equations 
in the regime of shallow water. With these expansions, we obtain a hierarchy 
of models for bi-layer shallow flows. First, we write lubrication models: using 
zeroth (resp. first order) expansion of the velocity field, we obtain a system 
of inviscid (resp. viscous) conservation laws on the fluid heights. We find 
the system of coupled Kuramoto-Sivashinsky equations in [TT] if we take into 
account of the capillary forces. We use this system of equations to study the 
spectral stability of steady states in the low frequency regime. Next, using 
first order expansions of the fluids velocities, we derive inviscid shallow water 
equations. For this latter step, one has to carry out a closure procedure: we 
have chosen to write the tangential stress at the bottom and at the interface 
proportionnal respectively to the average velocity at the bottom and the 
difference between average fluid velocities (up to correction terms). 



2.1 Description of bi-layer flows in the shallow water 
regime 

In this part, we write Navier-Stokes equations for bi-layer flows in a nondi- 
mensional form in the regime of shallow waters. We then perform an asymp- 
totic expansions of solutions with respect to the so-called aspect ratio (defined 
hereafter) in the neighbourhood of a Nusselt steady solution. 



2.1.1 Scaling Navier-Stokes equations 

Let us consider the superposition of two incompressible and immiscible fluids 
with density, viscosity and capillarity (pj, z/j, aj), z = 1,2 flowing down an 
inclined plane with a slope 9 (see flgure[2]). We introduce the aspect ratio e, 
the Reynolds number Re, the Froude number F and Weber numbers Wi as 

L z/i gH piHU^ 

where H denotes the characterictic depth of the fluid and L the characteristic 
length in the streamwise direction. The characteristic fluid velocity U can 



3 



(P2,i^2,(j2) 



hi{x,t) + h2{x,t) 



hi{x,t). 



X 



Figure 1: Two fluids flowing down an inclined plane. 



be chosen as the average velocity in the fluid layer for a Nusselt flow. We 
further introduce the additionnal numbers 

Pi ^1 

The motion of fluids (1) and (2) is described by Navier-Stokes equations 

Qi {dm + d^cUf + dmWi) + = ^ + {dzzUi + e^d:cxUr) , (1) 



F2 sRe 

Qi {dtWi + d^UiWi + d^wf) + -^-^ = + {dzzWi + e^d^^w^) , (2) 

dxUi + dzWi = Q, i = l,2. (3) 
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Here Qi — 1, g2 — p, l^i — 1, H2 — ^- These equations are set in the fluid 
domains 

^ ^ {{x,z) e <z<hi {x, t) } 

and 

^2,t = {{x,z) e R'^/hi{x,t) <z<hi + h2{x,t) = h{x,t)} . 

The kinematic conditions at the bottom, fluids interface and free surface are 
respectively 

ui{0,x) ^wi{0,x) ^0, ui{hi) ^ U2{hi), Wi{hi) ^ W2{hi) , 

(4) 

dth + ui{hi)d^hi = wi{hi), dth + U2{h)d^{h) = W2{h). 

We assume the continuity of the fluid stress at the fluids interface and at the 
free surface. First, the continuity of normal stresses yields 

P2{h) = r 7^d^U2ih) ^ ' 



Pi{hi) -P2(^i) = 



(i + £2(a,/i)2)i Re " ' ' i-£^{d,hy' 



-^{d^Mh)-r^d^U2{h))^-^,^^, 

with Ki — e^Wi. In order to take into account of the surface tension effects, 
we assume Ki — 0{1). Next, the continuity of tangential stresses yields 

{d:,U2 + e^d:,W2){h) = Ae^ — 9^U2{h) q 
v\d^U2^e d^W2){hi) - [d^ui + e d^wi)[h{) = As 1 _ g2(g y — ^^^i- 

Let us now describe the stationary solutions of this system. The velocity 
fleld does not depend on x and t. The fluid heights are constant hi{x, t) — h, 
h2{x,t) = 1 — h whereas the pressure is hydrostatic 

Pi{z) — c{h — z) + pc{l — h), yo < z < h, P2{z) — pc{l — z), 

and the fluid velocities have a parabolic proflle 

- - 

ui{z) ^ X{p{l- h)z + hz - yO<z<h, 

^ 

U2{z) = Xh{p{l-h) + ^) + ^[{l-h){z-h) - ^^^), yh<z<l, 



1 A • 1 n 1 A -R'^sm^ 

where A is a constant denned as A = — —5 — . 

In what follows, we will analyse the bi-layer flows in the neighbourhood of 
such steady solutions: this yields a natural scale for the characteristic fluid 
velocity U and thus the constant A has to satisfy an extra relation. If one 
choose the ratio between the total mass discharge rate and the total mass of 
the fluid then 



/ ui + p _ U2 = h + p{l - h), 

Jo J h _ _ 

h + p{l-h) 



A 



2 

t + 3ph\l -h) + Sp^hil -hf + ^{l- hf 



Note that for a single layer of fluid (p = z/ = 1), one recovers the condition 
of Vila A = 3. Another possible choice for the characteristic velocity would 
be the fluid velocity at the free surface: one then recovers the classical value 
A = 2. In both cases, there is a relation between the Reynolds and Froude 
numbers. Here, we have chosen U so that A = 3: there remains .Rg, 9^ h, k^, 
p, V as free parameters to design an experiment and describe a bi-layer flow 
of Newtonian fluids. 



We will derive shallow water equations from Navier Stokes equations inte- 
grated across each fluid layer (see e.g. [18], [16], [7]). First, we integrate the 
divergence free conditions on each fluid layer: using the kinematic conditions 
(jlj), we find the mass conservation laws: 

rh\ nh 

dth + dx{ ui{z)dz) = 0, dth2 + <9x( / U2{z)dz) = 0. 

Jo Jhi 

Denote qi = hiUi = J^^^ Ui and q2 = h2U2 = J^^ U2 the discharge rates in the 
streamwise direction: the mass conservation laws then read 

dthi + d^{hiui) = 0, dth2 + <9^(/i2M2) = 0. (5) 

Now, we write a system of equations which governs the evolution of = hiUi. 
This is done through the integration of momentum equations across each fluid 
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layer: 



eRe 



A . 2e ' 



dt{p / U2) +dj j 

J hi ^Jhi 



eRp Rf, 



hi + —dx / d^ui - T 







2 , P2\ , K2dxhdxxh 
PU2 + -F^)+ -J 



{i + e^idxhyy 

eRp 



2e ( f 
+ s-dx I dxi 

^ J hi 



U2]+T, (6) 



with T defined as 



In order to write this evolution system in a closed form, one has to find a 
relation between the different integrated quantities, the tangential stresses 
at the wall, at the fluids interface, T and the unknowns hi, qi. We follow the 
method introduced by Vila [T9] in the case of a single fluid layer. We expand 
the velocity field with respect to e in order to find an expansion of the above 
quantities and qi as functions of hi and their derivatives to any fixed order. 
For a given order, this enables us to write the unknown quantities in system 
^ as functions of (/ij, qi) and derive a shallow water model in a closed form. 

2.1.2 Asymptotic expansions of solutions to Navier-Stokes eqs 

In the shallow water regime e ~ 0, the fluid velocities and pressures Uj, Wi^pi 
almost satisfy a differential system in z. The "horizontal" fluid velocities 
Ui,i = 1,2 are solution to: 

eR 

PidzzUi + 0i>^ = eReQi{dtUi + Uid^Ui + widzUi) + -^d^Pi - PiS^d^xUi. (7) 
We add the boundary conditions: 

dzU2{h)=Ae'^ — "^^'^^/^^ 2 dxh - e'^dxW2{h), 
1 - e'^[dxh) 

n (^dxU2 — dxUi) (hi) 

vdzU2[hi) - dzUi{hi)=4e ^- dxhi 

1 - e^{dxhi)^ 

-e^ {udxW2{hi) - dxWi{hi)), 
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and Ui{0) — 0, Ui{hi) — U2{hi). 

The fluid pressures are solutions to the differential system 
dzPi+QiC = IJ'i^-dzzWi-e^F'^Qi[dtWi+UidxWi+WidzWij+iii^—dx^^ (8) 



whereas the boundary conditions for this system are given by 

P2[h) = 3- — d^U2{h)- , 

+ Re l-e'id,hy 

Ki F^d^^hi 



Pi{hi) -P2{hi) = 



Finally, the vertical velocities are solutions to 

dzWi = -d^Ui, wi(0) = 0, wi{hi) = W2{hi). 

There are three nondimensional numbers that are relevant to parametrize 
this set of equations: let us define 

« = P^eR^, S^'-0. (9) 

In what follows, we will assume that a, (3,6 ^ 1 so as to remain close to 
Nusselt type solutions. These assumptions are clearly satisfied when Re,F — 
0{1) but a wider ranger of parameters is valid. Next, we compute an Hilbert 
expansion of the fiuid velocity and pressure: 

oo oo 

k=0 k=0 

SO that 

k=0 k=0 

Let us first compute u[^\pf'^ Letting a,P,5 — >■ in the above equations 
leads a differential system in the z variable that is similar to the one which 
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determines stationary solutions. The fluid pressure is (up to this order) 
hydrostatic: 



pf\z) = c(p/i2 + hi- z) - KiF'^d^^hi - K2F'^dxxh, 

P2\z) = Pc{hl + h2- Z) - K2F'^droxh. 



The fluid velocities in the streamwise direction have a parabolic proflle 



u'i^ (z) = X {ph2Z + hiz- 

The computation of higher order terms is then straightforward: assume that 
we have computed uf\pf\ j < k, then uf''^^'' is calculated by computing 
the solution to 

Pid,,u^^'^^F,,k{u^^\p^^^), j<k, n = l,2 (12) 
with the boundary conditions 

dJ2'^'\h + h2) = gf \ udA'^'\hi) - d,uf^'\h,) = gt\ 



uf^'\h,) = u^t'\hi), ui'^'\Q) = 0. ^^^^ 
The solution uf^'^^'' to this system is 

uf"-'^ ^z{ugP - gP - F2,k{y)dy) - j^h,F,,,{y)dydz, 

u?^'^ = hi (u - - £ F2Mdy) - jT'' FiMdydz 

+ g'}\z -hi)- - [ I F2,k{y)dydz. 

^ J hi Jz 

We determine similarly an expansion of the fluid pressure to any flxed order. 



2.2 Lubrication theory 

The fluid velocities and pressures are expanded with respect to a, ^, 5 <S 1, hi 
and their space and time derivatives: we use zeroth and first order expansions 
of Ui to obtain respectively inviscid and viscous conservation laws for hi,h2- 
Then we study the spectral stability of steady states. 
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2.2.1 Inviscid and viscous conservation laws 



We first compute an inviscid system of conservation laws, which is the anal- 
ogous of Burgers equations in the case of a single fluid layer. The fluids ve- 
locities are given by Ui = uf^ + 0{a + f3 + 6). The discharge rates qi,i = 1,2 
are expanded as 



q2= l''u2 = \hh2{ph2 + ^) + — ^. + 0{a + (3 + 5). 



(14) 



2' u 3 

Inserting f|T4|) into the mass conservation laws (j5]) yields 
d,h + dJxhl{^ + ^))=0{a + ^ + 6), 
dth2 + (\h,h2{ph2 + y) + =0{a + P + 5). 

We drop 0{a + 13 + 5) terms in (|T5|) and obtain a system of partial differential 
equations for {hi, h2) in a closed form. A necessary condition of stability of 
steady state is that (|T5i) is a hyperbolic system. When this condition is 
satisfied, we obtain a useful information on the group velocities A_|_ > A_ 
of low frequency perturbations (respectively the velocities at the free surface 
and at the fluid interface): 

A± = ^(2p/i2 + Qpiyhih2 + Suhj ± Va) 
A = {hj (2p/i2 + hi) + Sp^hlhl) + ^phih\ (2p/i2 -hi)v + Ap^h^ 

Strict hyperbolicity is ensured if and only if A > 0. As it is a quadratic 
form in v the discriminant of A is — 128p^/if/i2 {hi + p/i2), which is negative 
when hi, /12, p are strictly positive. Then hi> Q and /i2 > ensure that the 
system (IT^ is strictly hyperbolic. 
Next, we use first order expansions of fluid velocities 



Uj = u. 



i^) + u^l) + 0{{a + P + 5)') 



to determine a more accurate system of equations. Inserting the expansion 
of qi,i = 1,2 into the mass conservation laws yields a system of Benney's 
equations (or Kuramoto-Sivashinsky when surface tension is considered) 



2^ ^ \hih2{ph2 + —) + ^hlj 



\(3d,(d{K)dJ^^ 



+ pd.(K{h)&'^(^^^),{lQ) 
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with the viscous coefficients dij defined as dij = — — dij^i — Xdij^2 

Re ' ' 

c^i,i,i = /ii(y + — ), c?i,2,i = p/^ily + y), 

phi , , 2 , ^? , P^2 7,2 , ^? 

"2,1,1 = ^ + P^1^2 + ^2ir, "2,2,1 = ^ + P"'l"'2 + P^2ir, 

h] /2h'\ 71pi/,o, 23p^z/,2,2 /P^ P^'^n, ,s P^ , 4 



phi /2u^hf 71pi/2 5pz/ 3p2z/2 Sp^z/, p ,, 

Z/2 V 15 120 ^ ^ ^ 24 4 ^ 1 2 -r ^ 1 2 -r ^ 2/. 

h2f5u^hl 25pz/2 4 Ilp2zy2 pj,^ 

— ;T /il ho + llt-llOr, 

1/2 V 24 24 1 2 ^ V g ^ g ; 1 2 



^2,1,2 = ^ ( ^ + ^/^?/^2 + (^ + ^-^)hlh^ 

5P'^.2.3 , ,V , 2pV,^ . , 2p3 



+ (p3.2 + ^),2,3^(^ + f|il),^,4^^,.)^ 
Ci2 2 2 = ^ f — + ^/.t/.2 + (^ + ^)/.?/.^ 

Z/2 V 24 ^ 24 1 ^ ^ ^ 2 2 ^ ^ ^ 

+ ^^hlhl + ' ' + -f^hl) . (17) 
The surface tension terms are given by 

n// \ ^1 ^2\ 7 2 /" ^1 ^2 \ 

^1,1 = /ii((«i + fi;2)y + '«2y), ^1,2 = «2/ii(y + yj, 

i^2,i = ^ + z.^/.^/^! + («:i + «:2)^, (18) 

K2h\ , , 2 K,2h2hl 

K2,2 = ^ + /t2/il/i2 + 

This system is in agreement with the one in [TTj. In the low frequency regime, 
the system f lT6|) of viscous conservation laws provides a criterion of spectral 
stability for the steady solutions which is consistent with the one given by 
Orr-Sommerfeld equations (if Re,F = 0(1)). One goal of this paper is 
to study the formation of roll-waves in bi-layer fiows. In the single layer 
case, they are the result ofthe competition between buoyancy and inertia. 
Therefore, we focus on the competition between inertia and buoyancy and 
their interaction with convective terms to describe the onset of roll-waves. 

2.2.2 Spectral Stability of Steady States 

Let us linearise f[T^ at a constant state {hi, /12): 

^* ( /I2 ) + '^^'^'^^^ {hl)= ^/^^(^O^- ( Jjj ) • (19) 
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Without loss of generality, we assume A/3 = 1. We have neglected the contri- 
bution of surface tension as they are not relevant in the low frequency regime. 
The dispersion relation is given by 

det{AId + ikJ{h) + k^d{h)) = 0, VA;gM. (20) 



and we assume \k\ <^ 1. We expand Aj,j = 1, 2 as Aj = ikAj. Equation (120|) 
then reads 

det(A + J(h)) = iktT(com{AId + J(h)fd(h)^ + 0{k'^). 

System f|T5|) is strictly hyperbolic : the eigenvalues Aj of J{h) are real and 
Ai > A2. Then Ai(fc),A2(A;) expand at = as 

_ iiicomUihi) -7i.jIdY d(hi)] 
AJk) = -ikA. = ^ + (21) 

As a result, stationary solutions are stable if 

tT(^com{j(h) -Aild)^d(h)^ < 0, tr (^com(j(75^) -A2ld)'^d(h)^ > 0. (22) 

We consider two particular situations: stable and unstable density stratifi- 
cation: we will see that in some particular situations, the Rayleigh Taylor 
instability may be suppressed by the convection. We also study the influence 
of inertia on stability properties. 

The spectral stability conditions fl22p have the simple form 

cotan^^ cotan^ . , , ■, 

ai{p)—B — <^a2(p), blip)— — >A62(p). 

ite lie 

It is easily seen that < 0: the free surface is then stable if 

Re < ^^f'^ cotan6'. 
Xa2{p) 

The situation is more involved for the fluid interface where bi can change 
sign. If bi{p)b2{p) < 0, the interface is stable when 61 (p) > and unstable 
otherwise. If 6i(p)62(p) > 0, the fluid interface is stable if 

bi{p){Re-^rr\cotane) <o. 

Ab2[Pj 
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In order to simplify the spectral analysis, we set /ij = 1, i = 1, 2 and consider 
the cases u < 1, u = 1 and u > 1 (stratification in viscosity). We have deter- 
mined stability curves Re = fk{p)cotaii9, k = 1,2 associated to the surface 
mode k = 1 and the interfacial mode k = 2. 

Case 1: < 1. We have chosen here u = 0.3, u = 0.7 and u = 0.9 that 
gives a good representation of all possible scenarii that arises as p varies. Let 
us first check the case u = 0.3. We have represented critical curves in picture 
[2j There exists pc ~ 3.3 above which 62 < otherwise both bi > 0,i = 1,2. 
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Figure 2: Values of bi,i = 1,2 and the critical curves fi(p) = ^— for 

^ ' J^\HJ cotan^ 

u = 0.3 

Then, for p > p^, the interfacial mode is stable whereas the full system is 
spectrally stable if Re < /i(p)cotan6'. If p < p^, fi{p) < fi^p)'- if Re is 
sufficiently small, the fiow is stable and as Re is increased, the surface mode 
is destabilized before the interfacial mode. 

Next we consider the case z/ = 0.7. The critical curves are represented in 
picture [31 There exists pc ~ 1.4 above which the interfacial mode is always 
unstable whereas the surface mode is stable for Re < /2(p)cotan^. If p < pc, 
there exists pi < p2 so that for any p < pi or p2 < p < pc, the scenario is 
identical to the previous case: as Re is increased, the surface mode is desta- 
bilized before the interface mode. If pi < p < p2, the interfacial mode is 
destabilized first as Re is increased. The case v = 0.9 is similar, except that 
for pi < p < P2, the interfacial mode is always unstable (see picture S]). 

Case 2: u = 1. There exists pc above which the interfacial mode is 
always stable. If p < pc, there exists pi < pc so that the interface mode is 
always unstable if pi < p < pc- If p < pi, the surface mode is destablized 
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Figure 3: Values of bi,i = 1,2 and the critical curves /i(p) 
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Figure 4: Values of bi,i = 1,2 and the critical curves /i(p) 
u = 0.9 



cotan6' 
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first Re is increased (see figure E]). 



Case 3: > 1. We have chosen u = 1.1 and u = 1.5, which gives 
a good representation of all possible scenarii that arises as p varies. We 
first consider u = 1.1 (see figure Ej). There exists pc ~ 3.5 above which 
the interfacial mode is always stable. Assume p < p^ and denote pi < p2 
respectively the first zeros of hi and the zero of 62- If Pi < P < P2, the 
interfacial mode is always unstable and if p < pi, the situation is similar to 
previous cases when p is small. If p > p2, the interfacial mode is stable only 
if Rp > /i(p)cotan6'. As a consequence, even at low Reynolds number, the 
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Figure 6: Values of bi,i 
V = 1.1 



1,2 and the critical curves /i(p) 
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interfacial mode is unstable. It is easily seen that there exists P2 < Ps < Pc 
such that the flow is always unstable if p2 < P < Ps, stable if ps < p < pa 
and 

/2(p)cotan6' < Re < fi{p)cota,n9. 

Finally, we consider the case u = 1.5 (see figure [7]). If pi is the zero of 
bl, then for any p > pi, the flow is always unstable under long wavelength 
perturbations. If p < pi, the situation is similar to the cases u < 1 when p is 
stable: the flow is stable at low Reynolds number and the flow destabilizes 
first through the surface mode or the interfacial mode. 
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Figure 7: Values oi hiA = 1,2 and the critical curves fAp) = ^— for 

^ •'^^^ cotan^ 

u = 1.5 

2.3 Shallow water theory 

The viscous conservation laws which govern the evolution of hi are sufficient 
to obtain a consistent stability criterion of constant states in the low fre- 
quency regime. However, the solutions to this system blow up in finite time 
when the fiow is unstable and may lead to some inaccuracy in the description 
of the motion of bi-layer fiows. In what follows, we consider shallow water 
models: indeed in the case of a single fiuid layer, they sustain nonlinear waves, 
so called "roll-waves" which are well known hydrodynamic instabilities. As 
a consequence, shallow water models are useful to describe the transition 
to instability in shallow fiows. Up to our knowledge, there's no consistent 
shallow water model which describes bi-layer fiows down a ramp (the single 
layer case was treated only recently ^16i|.^19j). 



Let us keep 0{1) and 0{e ^) terms in IQ: 

phi phi 

dt{ u,)+d,[j^ u\ 



[Xhi + vdzU2[hi) - dzUi[Q)) + 



Jhi ^ Jhi 



h 



— {Xph2-ud,U2{0)} —2 ■ (23) 

We first compute an expansion of the integrals: the integrals of the pressures 
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are given by 



/ Pi= pf^ + 0{a + l3 + 5), 2 = 1,2, 

Jhi-l Jo 

/ Pl^ = + PC/il/i2 - hiF'^ ((k2 + Kl)d.^xhl - K2dxxh2) , 

Jo, ^ 

/ = pC^ - K2F'^h2dxxh2- 

whereas the integrals of convection terms are given by 

u\= \ (nfV + 0(a + /3 + 5), ^ = l,2 

1 Jhi-\ 



(24) 



2 



+AH;(A5 + pft^fti + ^)). 

In both cases, we only keep 0(1) terms. In order to write these quantities 
as classical convective terms, we introduce Qi^i = 1,2 so that 

C\uff = hul + Qi(/ii, h2), j\uff = h2ul + Q2{hl, h2). (25) 

Jo J hi 

Here, Qi depend on hi and are defined as 

.2P^ ^2 



Qi = A^/^?(| + ^(/.i+pM), g2 = A^ 

Inserting and (12^ into (1251) . one finds 

hi 



45 1/2 



/ C ft \ 

dt{hiUi) + (9^(^/iiMi + —{phih2 + y) + j - hidxxx{i^2{hi + /la) + ^i/ii) 



= ^h2dxhi + 1^ + -^{i'd:,U2{hi) - 9^Mi(0)), 

p[dt{h2U2) + dx{h2ul + + Q2) j - i^2h2dxxxh 

= -^h2dxh + - -^d,U2{hi). 

This system is almost in a closed form. Let us now write dzUi{0) and 
dzU2{hi) as functions of Ui and hi. We clearly see that an expansion of Ui up 
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to order 0[{a + (3 + 5)^) is needed. We use the method introduced by Vila 
[TU] in the case of a single fluid layer to write these terms in a closed form. 
In [19], the wall stress is chosen proportionnal to the average velocity. One 
has to expand both the wall stress and the average velocity up to order 1 
and obtains an expansion of the wall stress with a zeroth order term propor- 
tionnal to the average velocity, the next term depending on the fluid height 
and its time and spatial derivatives. See [1] for a mathematical justification 
of this derivation. 

The situation is more involved for bi-layer flows and several closures are pos- 
sible. However, in order to fit with the model in [19], we search an expansion 
of the fluid stress at the bottom in the form 

dzUi{0) = 71(^1,^2)7^ + (h.o.t), 
hi 

with 7i(/ii, 0) = 3. The fluid stress at the bottom is given by 

dM0) = liihi,h2)^ + Ri, (26) 
til 

with 7i and Ri defined as 71 = 6— -^--^—^—7— and 

2hi + 3pn.2 

^1 = ~ nt. ^,^0 I. {Ri,id^h + Ri,2dxh2) - 5c{p^ - p) ^J^^^"^ , 19^/^2 
2hi + 3p/i2 2hi + ph2 

+ oL^^^o\ ^xxxif^iphi + K2{p - l)ihi + /la)), 
Zhi + Spn2 

^ = —h\ + ^h\h2 + ^-^hlhl + (^ + ^)hlhl + ^hihl 

' 15 ^ 60 ^ 10 ^ 2 ^ 4 3z/^ ^ ^ 3z/ 

,= '^hl + + (^ + ^)hlhl + ^hlhl + -^hihi 

^'^ 15 ^ 60 ^ ^ ^ 20 4i/^ ^ 2 12z/ ^ 2 3z/2 ^ 



Next, we write the fluid stress at the interface as 

U2 - Uint 



dzU2{hi) 



ho 



with Uint = ui{hi) = U2{hi) and expand Uint as Um = l2{hi, h2)ui + (h.o.t). 
The fluid stress at the interface reads 

3 

dzU2{hi) = —{U2 - 72 Ml) + R2, 

h2 
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with 72 and R2 defined as 72 = ^—r~~~^~r~ 

2hi + 3p/i2 



= -^/.? + ^hlh2 + ^^h\hl + + ^-^)hlhl 



20 ' 5 ' 10 ' ^ '2i/ 4 
+ (#^ + ^)/i?/i^ + :^/ii/ii + 



20 ' 5 ' '8i/ 120 ^ 4i/ 
■^2hih2 + {^2+^)h^h, + —l 



Note that we have imphcitely used the mass conservation law 

dthi = + (h.o.t) 

to transform time derivatives into spatial derivatives. As a result, we obtain 
a shallow water model for bi-layer flows in a closed form: 

dthi + d.,{hiui) = 0, dth2 + dx{h2U2) = 0, (27) 

dt{h,u,) + dx{h,ul + ;^) = ^(A/M + 3^^^^^ - 371^) + nm 

Ir ^ £ Re ri2 111 

d,{h2U2) + dx{h2ul + ^) = ^[Xh2- ^-^"^IZJ^) + (29) 

Zf £ rig p fl2 

c ~ c ~ 

T^i ^ -^h2dj;hi + TZi, 7^.2 = --^/ii5x^2 + ^2, 

where TZi, i = 1,2 are only fonctions of hi and their spatial derivatives. 
These are corrective terms to the hydrostatic repartition of pressure within 
the fluids which are due to surface tension, buoyancy and inertia. They are 
written as 

TZi = -d^Qi + ^(i/i?2 - Ri) + hidxxx{K2{hi + /i2) + ki hi), 

£Ke 

T^2 = — dxQ2 ^-^2 + —h2dxxx(hi + /i2). 

p £ReP p 
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This system is not in a conservative from and this may lead to some inde- 
termination in the presence of shocks. One can drop this indetermination by 
the use of "nonconservative paths" (see the next section for definitions). 
For convenience we rewrite the system in matrix form. Set W = {hi, h2, ui, U2Y' ■ 
We write system (l27l[28| [29!) as 



dtW + d^FiW) = B{W)dxW + G{W), 
where F is defined as 



(30) 



F{W) 



hi 



chl 
2F2 



whereas B is given by 



+ \'hl 



Qi 

Q2 

^+'-^{hi+ph2 

45 12 ^ ' 



Q2 chl 
h2 2F2 



+ 



45z/2 



B{W) 














" 




dxhi 


dW 
















dxh2 


dx 


Bii 


B12 










dxQi 




B21 


B22 










dxq2 



with 



Bii 



-^h2 + 



F2 2hi + 3ph2 



h2 {2hi + 3ph2 



7^2,1 



B 



12 



h2 
A2 



3/^2 

R2,2 + C{p'- p) ^2 



{2hi + 3ph2) 
hih2 



^2hi + 3ph2^''' + T^' " 2hi + ph2 



B. 
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vX^ 



h2p {2hi + 3p/i2) 



-^2,1 



B' 
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V 



p{2hi + 3ph2) F2 



hi. 



The source term G (W) reads 
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G(W) 



1 



Xhi + Su 






U2 - I2U1 




h2 

6u U2 - 72M1 



A/12 + 



P h2 



3 Roll-waves in shallow water equations 

In this section, we prove the existence of roll-waves in bilayer flows when 
steady states are unstable. They are defined as piecewise smooth and spa- 
tially periodic travelling waves, entropic solutions to shallow water equa- 
tions. For hyperbolic conservation laws, the shocks must satisfy the Rankine- 
Hugoniot and Lax shock conditions. In [6] , these solutions are proved to exist 
in a single layer of fluid modeled by a shallow water system. For general hy- 
perbolic conservations laws with source terms, there are small amplitude 
roll- waves [14j. We generalize this result to nonconservative hyperbolic sys- 
tems in order to deal with our shallow water model for bi-layer flows. We 
also carry out direct numerical simulations to show the existence of large 
amplitude roll-waves. 

3.1 Existence of small amplitude roll- waves 

In this section, we consider the problem 



We assume that system (pTj) is strictly hyperbolic in the neighbourhood V(mo) 
of a constant solution u = uq, meaning that A{u) has n real distinct eigen- 
values (Afc(M))^^^ ^, 

Xi{u) < ... < Xk{u) < ... < A„(n), Vm G V(uo). 

We suppose that both A and g have a power serie expansion at u = uq with 
a disk of convergence containing B{uQ,r) for a suitable r > 0. As in [6], 
we search a spatially periodic travelling wave u{x,t) = U{x — ct) with U a 
2L-periodic function with discontinuities at Xj = (2j + j & 1^ which 
satisfies the differential system 



du 
'dt 



du 
dx 



x e M, t>0. 



(31) 



{AiU{x)) -c)U' = giU{x)), Vx e L). 



(32) 
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Next, we formulate conditions for shocks at Xj = {2j + 1)L, j G Z. For that 
purpose, we need to define the so-called "family of paths" [5j. These paths 
were introduced to give a rigorous definition of nonconservative products in 
hyperbolic systems. A family of paths $ in C is a locally Lipschitz 
map $ : [0, 1] x x — )■ fi, such that 

• $(0,M;,Mr) = Ui and (t>{l,Ui,Ur) = Ur, foY any Ui,Ur G 

• for any bounded subset O G Q, there exists a constant k such that 
-^{s,Ui,Ur) < /c I Mz — Ur- 1, for all tij,, G O and almost every s G [0,1]. 

for any bounded subset O G Q, there exists a constant K such that 



(9$ 
ds 



9$ 

ds 



< K{ 



\u, 



U, + 



for all G O and almost every s G [0, 1]. 

When such a family has been chosen, one can define generalized Rankine 
Hugoniot jump condition across a discontinuity with speed ^ 

^ {^ild - A{^{s,u~ ,u^))^^{s,u~ ,u+) ds = 

where u^,u^ are the left and right limits at the discontinuity. For the prob- 
lem of roll-waves, this can be written as a nonlinear boundary condition, 
setting = U{L) and u+ = U{—L): 

(^cld- A{<^{s,U{L),u{-L)))^^{s,U{L),U{-L))ds = 0. (33) 

This Rankine- Hugoniot condition is completed with a Lax shock condition 

Xk{U{-L)) < c < Xk{U{L)), Xk-i{U{L)) < c < Xk+i{U{-L)), (34) 

for some k, 1 < k < n. Herein, a roll-wave is a solution of the so called 
"roll- wave problem" fl32ll33f34p . 

Let us fix /e, 1 < A; < n: we denote rk{u) the eigenvector of A{u) associated 
to the eigenvalue Xk{u) and we assume that the characteristic field rk{u) is 
genuinely nonlinear VXk{u).ri:{u) ^ 0, Vm G V(mo). We define Hfc('u) as the 
projection onto Ker(^(M) — Xk{u)Id) with respect to lm{A{u) — Xk{u)Id). 
The eigenvalue is isolated so that u i— )■ Iik{u) is C^. Ker(^(M) — Xk{u)Id) is 
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one dimensional and we identify, for any v G M", Ilkiujv to a real number. 
Finally, we will suppose that $ is analytic in each variables in order to 
simplify the discussion: this hypothesis is clearly statisfied for straight lines, 
a natural choice in numerical schemes. Let us prove the existence of small 
amplitude roll- waves. 

Theorem 1 Assume that 

Uk{uo)DA{uo).rkiuo).rk(uo) 7^ 0, Uk(uo)dg(uo).rkiuo) ^ 0, (35) 

T-r X n /i^- ^ 7=^\ 7^^^k{uo).rk{uo) > 0, 36 

Uk{uo)DA{uo).rk{uo).rk{uo) 

and dgiuo) : M" — > M" is invertible. Then there exists a family of small 
amplitude roll-waves solutions of I^SEM^y parametrized by wavelength. 



Proof. We will prove that the solutions to fl32|33|34p are zeros of a submer- 
sion between suitable functional spaces. First, let us recall the construction 
of a formal roll-wave. Set L = 2riT. We search a roll-wave in the form 
u{x,t) =u + r]v{^), 1] <^1. The system fl32ll33ll34p reads 



[A{u + r]v{x)) — c)v'{x) = T g(u + r]v{x)), Va; G (— 1, 1), (37) 

J ^c/c?-^($(s,M+77t;(l),M+r7i;(-l))) j^(s,II+r7t;(l),M+r7t;(-l)) ds = 0. 

(38) 

As 7] —)■ 0, the Lax shock conditions are 

Xk{u + < c < Xk{u + vv{^)). (39) 

Letting r/ ^ in ([38]) yields A{u){v{l) - v{-l)) = c{v{l) - v{-l)). Nec- 
essarily c = Xk{u) and v{l) — f (— 1) is an eigenvector of A(u) associated to 
\k(u). Let us search v{x) = a{x)rk(u): the Rankine Hugoniot condition is 
then satisfied. If u = uq a zero of g, then dividing f p7|) by rj and letting 
— > yields 

DA(uo).rk{uo).rk(uo) a{x)a'{x) = dg(uo).rk(uo)a{x). (40) 

Then, a projection onto Ker(^(Mo) — Xk{uo)) yields a{x) = Tx. As 77 — )■ 0, 
the Lax shock conditions are 

VAfe(Mo)-rfc(Mo)a(-l) < < V \k(uo)-rk(uoMl)- (41) 
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This relation is clearly satisfied under the assumption fl5^ and the construc- 
tion of a formal roll-wave is complete. In order to prove the existence of 
roll-waves close to this formal solution, we introduce 

j;,^ : X X V(mo) Y X R" 



+ (Jrf - nfc(n)) ((^(U + r^t;) - Afc(M))t;') 

+ r(/(i - nfc(M))( < giu + r]v) > -g{u + r]v)), 

'''^ A{u) - A{u + riv{x)) 



J-'ri,T{v,u)2= / g{u + riv{x))dx + T] / v'{x)dx 

J~i J-1 V 

+ A{(l>{s,u + r]v{±l)) - A{u))^{s,u + T]v{±l))ds, 

with < giu + T]) >= |J%;,r(w, f )2- The functional spaces X, Y are defined as 

Xo = I / G C\-l, l)/fix) = ^ a„ x", ^(n + 1) 

L n>0 n>0 

/eXo/(i-nfc(Mo)) j fix)dx = o 



aJ < oo 



X 
Y 



I / e C\-l, 1)/ fix) = ^ a^x", ^ |a„| < oo I 

L n>0 n>0 J 



The operator J'ri,T is well defined and C^. It is clear that a zero (m, v) of J-"^,,- 
corresponds to a roll-wave (see [H] for more details). Let us fix tq > 0. As 
— > 0, it is easily seen that a zero (u, v) of J^o,ro satisfies 

1 /•! 

nfc(n)Z}^(M).t>(x).f'(x) — nfc(u)(i(7(M) (t>(x) — - v{s)ds), 

(42) 



(1 - Uk{u)) (^A{u) - Xk{u))v'{x) = 0, Vx e (-1, 1), 
g{u) = 0. 

The roll-wave (uq, v{x) = a{x)rf:(uo)) is clearly a zero of J^o.ro- The end of the 
proof is similar to the conservative case. One prove that DJ^o^roi^o, C({x)rk(uo)) 
is invertible and the implicit function theorem applies: for < ?7 ^ 1 and 
T ^ tq, there exist a unique zero of J^r],T which is close to the formal roll- wave. 
If condition f l36p is satisfied, the roll-wave satisfies Lax shock condition 
for ?7 > sufficiently small and the proof is complete. □ 
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A slight modification of tliis argument enables us to deal with "physical" 
source terms (which is the case for bi-layer flows). Suppose that g has the 
particular form g{u) = (Oirp,/i(m)) and dh{uo) : M" — > M"^^ is onto. The 
construction of the formal roll-wave is the same. Then one can prove that 
for < ?7 -C 1 and r ~ tq, there exists a family of roll- waves solutions that 
belongs to a p + 1-dimensional manifold. In this case, J^r],T is a submersion 
at the point corresponding to the formal roll-wave. 



3.2 Numerical simulations 

In this section, we investigate numerically the existence of roll- waves through 
direct numerical simulations of the shallow water equations fl27 t [2^. We 
consider the case where all the interfaces are unstable (these latter solutions 
are the bilayer counterpart of regular roll- waves into a single fluid layer). 

3.2.1 Numerical scheme 



We use a classical upwind difference scheme as described in [13] . We assume 
X G [0, L] and integrate (pi on the time interval [0, T]. System (15U]) is strictly 
hyperbolic provided that the eigenvalues of M [W) = A {W) — B {W) are real 
and distincts. Here A {W) denotes the Jacobian matrix of F: 



A 



OF 
dW 









1 














1 


^31 


^32 


hi 








Ai2 





2g2 
ho 



with (recall that the nondimensional numbers /9, 5, A are defined as /3 = eRe 
^^ = ji and A = -^sin^) : 



^1 

^32 
^42 



X'hl 



P_ 
12 

ch2 



{hi + ph2] 



12 



hi F2 



5AV^ 
45z/2 



p\^h\ 



h 

12 



+ \'hl 

ph2 
6 



2/^1 ph 
45 12 



As a consequence the matrix M is given by 
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M{W) = A{W) - B{W) 









1 














1 


M31 


M32 


2qi 
hi 





M41 


M42 





2g2 

/i9 



and 



M- 



31 



2/^1 p/i2' 
45 12 



/l2 (2/li + 3p/i2 



32 



+ 



/ 2 \ 3/1? 
2/ii + 3p/i2 F^""^^ ^^2/ii + p/i2 



A 



hi ph2 

12 6 

2 



41 



42 



--^2,1 



/12P (2/ii + 3p/i2) 

^2 c/l2 5A^p/l2 C 



/l2 + ^2 + 45j,2 



-R2,2 + C (p' - p) 



3^ 

2F2 



p (2/ii + 3p/i2 



We approximate system fl5Ul) by the regularized system: 



SW" 9 dW 
— ^-F(W)^G(W) + BiW)^ 



dx 



dx 



(44) 



Ax 

where (W^) •^ajVT) represents a numerical diffusion introduced by 

the scheme. The diffusion matrix D (W) must take into account that the 
effective transport term in (150]) is M {W) d^W, where M is given by We 
propose to discretize the nonconservative product using the relation BWx = 
(BW)^ — BxW. We use a two-stage second order time scheme. 



26 



The numerical scheme is then written as 



fiwn=Giw.p 



Ax 



2Ax 
Ax 



where is the second order MUSCL reconstructed state of using clas- 

sical flux hmiter function minmod as described in pD], = — - is 

an intermediate state between VT" and W^i, and the numerical flux 0" is 
given by 



Fc {W-, W-^,) - -^D (W-, WJX, 



Ax 



where Fq and D are respectively approximations of F and V ai x = x^.i 



2 



D{U,V)=X\A\X-^ 



with A the matrix of the eigenvalues of M (^^-j^) and X the matrix deflned 
by its eigenvectors. Note that we note (c?i)j=i 4 the eigenvalues of D, the 
CFL condition is then given by : 

At 

- — max di < 1 

Ax i=l,...,4 

3.2.2 Numerical results 

We carry out numerical simulations when the steady states are unstable. We 
have flxed: z/ = 0, 9 and p = 0,5 (see flgure H]) and set 6* = 7r/4 so that 
Re = f and = -\/2//6. The aspect ratio e is 0, 01 <^ 1. We made numeri- 
cal simulations for / = 0, 5 (all interfaces are unstable). The initial condition 
is the steady state W = (1, 1, 1.75, 3.56)"^ perturbed in the direction of the 
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most unstable eigenvector. 



Recall that we fix / = 0.5 and (-Re = /, F"^ 



V2f 



6 



). Let us build the 



initial condition: the unique unstable eigenvalue of matrix d^S — 2ttM is 
A2 = 13.4444 + 3.075? and the associated eigenvector is 



A, 



0.0436415 - 0.0321769Z 
-0.0916431 - 0.400677Z 
0.047492 + 0.109129? 
0.902194 + 0? 



We start the numerical simulations with the initial condition 

Winit = W + 5.10-3 (cos {2ttx) - sin (27rx) $7) . 

We took 250 points for one period in the spatial mesh. At t = 0, the fluid 
interface and the free surface are periodic with the same period but different 
amplitude. As we used a 50/00 perturbation, interfaces are close to steady 
state. Note that the scale is different for internal and free surface wave: on 
the left is the scale for hi + /i2 whereas the scale on the right corresponds 
to the free surface hi. We clearly see the formation of roll- waves both at 
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Figure 8: Fluid heights at time t = 10 
the free surface and at the interface and that they are in phase. We have 
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computed the spatial Fourier transform of this signah it is composed of 60 
different modes and the first 20 modes are the most relevant. We also plotted 
in picture [9] the time evolution of the first two Fourier modes : for t G [0, 1.2], 
the amplitudes of both interface does not vary much and for t G [1.2, 2], there 
is creation of roll-waves which then stabilize. 



Time evolution of amplitude of the two first modes 
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Figure 9: Time evolution of the first two Fourier modes 



4 Conclusion 

In this paper, we have obtained consistent shallow water equations for bi- 
layer fiows from the Navier Stokes equations in the presence of capillarity. As 
a byproduct, we carry out a complete spectral stability analysis of bi-layer 
fiows in the low frequency regime. We proved that this system is a gener- 
alization of the system of Kuramoto Sivashinsky equations derived in [TT] 
and that it is useful to describe nonlinear waves of arbitrary amplitude in bi 
layer fiows. Numerical simulations then confirm the existence of well known 
hydrodynamic instabilities, so called roll-waves, which could be localized on 
the fiuid interface or on both interfaces. 

This system of shallow water equations is a hyperbolic system in a non conser- 
vative form, a common property in shallow water systems describing bi-layer 
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flows. Therefore, tliere is non uniqueness in tfie definition of shocks. One 
possibihty would be to derive higher order shallow water models with a van- 
ishing viscosity: the physical viscous term would then select the "physical" 
jump conditions. For applications purposes, it would be also of interest to 
derive bi-layer models for non Newtonian fluids in the spirit of [3] , where con- 
sistant shallow water equations for single thin layers of Bingham and power 
law fluids. 
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